scm6A-seq reveals single-cell landscapes of the dynamic m6A during oocyte maturation and early embryonic development

N6-methyladenosine (m6A) has been demonstrated to regulate RNA metabolism and various biological processes, including gametogenesis and embryogenesis. However, the landscape and function of m6A at single cell resolution have not been extensively studied in mammalian oocytes or during pre-implantation. In this study, we developed a single-cell m6A sequencing (scm6A-seq) method to simultaneously profile the m6A methylome and transcriptome in single oocytes/blastomeres of cleavage-stage embryos. We found that m6A deficiency leads to aberrant RNA clearance and consequent low quality of Mettl3Gdf9 conditional knockout (cKO) oocytes. We further revealed that m6A regulates the translation and stability of modified RNAs in metaphase II (MII) oocytes and during oocyte-to-embryo transition, respectively. Moreover, we observed m6A-dependent asymmetries in the epi-transcriptome between the blastomeres of two-cell embryo. scm6A-seq thus allows in-depth investigation into m6A characteristics and functions, and the findings provide invaluable single-cell resolution resources for delineating the underlying mechanism for gametogenesis and early embryonic development.

experiments would be necessary to support the authors' conclusions.
Major comments: 1. Fig. 1D shows that the m6A modified fragments detected by scm6A-seq are the highest enrichment ratio in CDS region, which was not consistent with the previous reports that m6A was mainly found in 3'UTR, as mentioned in the Introduction (line 53). Please explain this difference. Supplementary Fig. 1e, more than 10000 m6A modified RNAs are detected in the groups of 50 and 70 GVs by bulk MeRIP-seq, while, the mounts of modified RNAs are about 5000 in the group of 25 GVs by the same way and so it is with single-cell through scm6A-seq. How to achieve the optimal balance between cell number and detection amount? 3. Some statements lack for sufficient evidence in the article, for example, Supplementary Fig. 1h suggests that scm6A-seq is more sensitive than MeRIP-seq through the expression and modification of Hist4h4 gene; the integrated genomics viewer tracks of Lars2 and Mtus1 show the commonly and specifically m6A modified RNAs at different oocyte stages in Supplementary Fig. 4e, but it does not reflect the "far more" described in line 288. Please provide more evidence to support these descriptions. 4. At present, some studies have detected m6A modified RNAs at the developmental stages of oocyte and pre-implantation embryos, such as references 6 and 10 cited in this paper. Has the scm6A-seq data gained from this study been compared with these publicly available data? 5. As observed in the study, m6A modification plays the opposite regulatory role at different stages, promoting mRNA degradation in GV oocytes while protecting the modified RNAs from degradation during MII-to-zygote transition. This phenomenon is interesting, please explain it further. 6. I would like the authors to reconsider the Discussion. For example, some existing studies on RNA m6A modification in the field of reproduction should be included in the Introduction section in Introduction. In addition, most of Discussion are merely repetition or summary of the results.

As shown in
Minor comments: 1. The ordinate of Fig.1F is marked incorrectly.
2. Line 157-158 has a syntax problem. 3. Fig. 2e does not reflect the role of mettl3, although it is described in the legend, please describe it as consistently as Fig. 4e Revision Summary

Summary of major comments from Editors and Reviewers.
Based on the requests and comments from the reviewers, we have performed further analyses and designed additional experiments that have largely improved our manuscript. The revised parts have been marked with red color in our revised manuscript and figures. The detailed revision results are summarized and listed in the following Table 1.  Fig. 8a, Revised Fig. 2e) using public data 5 .
Combined with the public eCLIP-seq data 6 , we also found the METTL3-dependent m 6 A modified peaks are close to the YTHDF2 targets (Revised Fig. 2e) and more stable in Mettl3 Gdf9 cKO oocytes (Revised Fig. 2c). These results suggeste that YTHDF2 is the probable regulator for the m 6 A mediated RNA degradation in the GV oocytes.
For the relative higher stability of m 6

5
The statements about the sensitivity of scm 6 A-seq and the stage specific m 6 A methylome lack sufficient evidence in the article.

2#: Q3
To further illustrate the sensitivity of scm 6 A-seq for m 6 A detection, we have compared scm 6 A-seq and bulk MeRIPseq for GV oocytes and found that nearly 70% methylated RNAs in 2500 GVs bulk MeRIP-seq were also identified in scm 6 A-seq (Rebuttal Fig. 6a). Besides, there were 50% methylated RNAs only detected in single GV by scm 6

6
The comparison of the scm 6 A-seq data gained from this study with the publicly available data.

2#: Q4
To compare the scm 6 A-seq data with previous relevant publications, we have analyzed bulk MeRIP-seq data for GV oocytes with two replicates 8 , and low input ULI-MeRIP-seq data for GV oocytes with three replicates in another recent publication 7 . We found that the modified RNAs detected by scm 6 Fig. 1a, Revised Fig. 2d).

Point-by-point responses to reviewers
For the METTL3-dependent m 6 A modified RNAs, we have elucidated that these RNAs tend to be degraded in GV oocytes (Revised Fig. 2g; Revised Supplementary Fig. 2e,f), thereby showing globally increased expression in Mettl3 Gdf9 cKO oocytes. For the METTL3-independent m 6 A modified RNAs, the globally increased expression may be due to the association of m 6 A with RNA abundance. Thus, we have further evaluated the abundance of both m 6 A modified and unmodified RNAs in wild-type and Mettl3 Gdf9 cKO oocytes, and found that the m 6 A modified RNAs display a higher expression level than the unmodified ones (Rebuttal Fig. 1b, Revised Supplementary Fig.  2d), which is consistent with the previous study performed in single cells using scDART-seq 1 . In addition, we have also analyzed the data from a previous study on KIAA1429 2 , one known component of the m 6 A methyltransferase complex, and found that the regulation of m 6 A in promoting RNA degradation in GV oocytes was also observed in Kiaa1429 cKO oocytes (Rebuttal Fig. 1c).
Collectively, these data demonstrated that m 6 A modification contributes to RNA degradation in mouse oocytes. According to the current reports, the m 6 A modification could be catalyzed by different methyltransferases, including METTL3 and METTL16 9 .

The conclusion in Fig5 that m 6 A protects its modified RNAs from degradation during MII to zygote transition is not such strong in Fig5b and 5c. Can you use the RNA-seq data of Mettl3 cKO oocytes which will arrest in zygote to supplement the data?
Response: Thanks for this thoughtful suggestion. To comprehensively analyze the regulation of m 6 A on RNA metabolism during MII to zygote, we have first analyzed the expression changes of the unmodified RNAs and found that compared to the increased expression level of m 6 A modified RNAs (Original Fig. 5b, Revised Fig. 5b), the unmodified RNAs showed a significantly decreased expression from MII to zygote (Rebuttal Fig. 2a, Revised Supplementary Fig. 5f). Furthermore, we have followed the reviewer's suggestion and collected the Mettl3 Gdf9 cKO MII oocytes and zygotes for total RNA-seq. Strikingly, we found that the m 6 A modified RNAs (Original Fig. 5b,  Revised Fig. 5b) degrade more rapidly compared to the unmodified ones in Mettl3 deficient MII oocytes and zygotes (Rebuttal Fig. 2b, Revised Fig. 5d). We have also included this information in the Revised manuscript (Page 10, lines 279-291).
Rebuttal Fig. 2. m 6 A protects m 6 A modified RNAs from degradation during the MII-to-zygote transition.
a. Box and whisker plot showing the expression of unmodified RNAs in oocytes from MII oocyte to 4-cell stage. P values were determined by Wilcoxon test. b. Cumulative fraction of expression changes of m 6 A modified RNAs (blue) and unmodified (grey) RNAs in Mettl3 Gdf9 cKO MII oocytes and zygotes. P value was determined by Wilcoxon test. *, P < 0.05; **, P < 0.01; ***, P < 0.001.
3. The content of the whole manuscript are too much and complicated. Some information about the expression profile such like line 249-264, 285-296 may be removed to highlight the major points.
The subhead of each paragraph should be revised to strengthen the conclusion. Response: Thanks for these suggestions. We have followed the reviewer's advice to remove the content of expression profile of early embryo development and to highlight the major points in the Revised manuscript (Pages 9-10, lines 267-291). Additionally, we have reorganized the subheadings to strengthen the conclusion for each part in the

Line 145-146, add a brief description of the phenotype of Mettl3 cKO oocytes.
Response: Thanks for pointing out this. As "In Mettl3 flox/flox ; Gdf9-cre (Mettl3 conditional knockout (cKO)) mice, the size of the ovary largely decreases, and the oocytes are arrested from maturation. Few MII oocytes can be observed, and most of the oocytes are arrested at the GVBD stage. Moreover, the Mettl3 Gdf9 cKO oocyte-derived zygotes are arrested at 1-cell stage after fertilization". We have added this description about the Mettl3 Gdf9 cKO oocytes in the Revised manuscript (Page 6, lines 151-155).

In sup Fig3g, the control group are the total RNAs, it's better to choose random RNAs with the same gene number as development-silencing RNAs.
Response: Thanks for pointing out this. We have randomly selected the same number (n=1478) of development-silencing RNAs as a control group to further evaluate the expression change. We still observed a significant difference between randomly selected RNAs and development-silencing RNAs, but not between the total RNAs and randomly chosen RNAs (Rebuttal Fig. 3, Revised  Supplementary Fig. 3g).
Violin plot showing the expression change of total RNAs and the development-silencing RNAs between Mettl3 Gdf9 cKO oocytes and control oocytes. The same number of the developmentsilencing RNAs was randomly selected from the total RNAs as non-sense control. P value was calculated by unpaired t-test. 4. The order of Fig4e and Fig4g seems to be wrong according to the text. Response: Thanks for pointing out this. We apologize for this error and have corrected it in the Revised manuscript (Page 30, lines 921-922, 925-927) and Revised Fig. 4.

In Fig5b, how many genes are plotted? Number should be showed also in Sup Fig5g.
Response: Thanks for pointing out this. The detailed gene list and expression information have been included in Original Table S7. We have followed the reviewer's suggestion and added the numbers of modified and unmodified RNAs in each stage into the corresponding Revised Fig. 5b-c and Revised Supplementary Fig. 5g.

In Sup Fig6a, the variance among cells at the same stage seems too large to get the conclusion that the m 6 A level are decreasing after zygote.
Response: Thanks for the comments. We have modified this statement to "We observed a trend of decreased m 6 A level after ZGA, which is consistent with the result detected by HPLC 7 " in Revised manuscripts (Page 10, lines 305-306). Moreover, the variance may be caused by the highly heterogeneous blastomeres of early embryos. Our result showed that the ZGA process is asynchronous in the 2-cell embryos, and the variance of m 6 A level decreases afterwards (Original Fig. 6f,g). Thus, we draw the conclusion of "m 6 A decreases after ZGA".

Reviewer #2 (Remarks to the Author): In this study, Yao and colleagues develop a new technique of single-cell m 6 A sequencing (scm 6 Aseq) to profile the m 6 A methylome and transcriptome simultaneously. The authors observe the dynamic changes of m 6 A modified RNAs during oocyte maturation and early embryonic development. They find that m 6 A deficiency led to aberrant RNA clearance and low quality of oocytes in the model of Mettl3 conditional knockout mouse. The asymmetric epi-transcriptome that is dependent on m 6 A further explains the heterogeneity between the blastomeres of 2-cell embryo. Overall, this study is carried out by the method of scm 6 A-seq in the single cell dimension to investigate in-depth into m 6 A characteristics and functions at the stage of gametogenesis and early embryonic development. Comments for the author However, before publication, further interpretation and discussion of the data as well as additional experiments would be necessary to support the authors' conclusions.
Major comments: 1. Fig. 1D shows that the m 6 A modified fragments detected by scm 6 A-seq are the highest enrichment ratio in CDS region, which was not consistent with the previous reports that m 6 A was mainly found in 3'UTR, as mentioned in the Introduction (line 53). Please explain this difference.
Response: Thanks for pointing out this. To clearly explain the difference of m 6 A distribution, we first have searched the literatures about m 6 A distribution in human cancer cell lines, and found that this variation is probably derived from different methods. For m 6 A-SAC-seq, more than half of the sites are located in the 3' UTR region, and over 40% of m 6 A sites are located in the CDS region of polyA+ transcripts 10 . For DART-seq, there are nearly 70% of m 6 A sites in the 3' UTR region and about 20% in the CDS region 11 . For m 6 A-REF-seq, almost 70% of detected m 6 A sites are located in CDS, and about 30% of the m 6 A sites are located in 3' UTR 12 . Notably, for antibody-based methods, the distribution of detected m 6 A sites or peaks can be slightly influenced by different commercial sources of antibody 13 .
Another important factor affecting the distribution of m 6 A is the type of RNA for detection. There are much more m 6 A sites distributed in promoter, intron, and intergenic regions in total RNA samples than the polyA+ RNA samples 14 . Therefore, m 6 A can exist in any location of the transcriptome and is more enriched in long coding exons and 3' UTR of polyA + RNAs. We have revised the description in the Introduction section in the Revised manuscript (Page 3, lines 52-53).
Additionally, according to the previously reported methods 8,15 , we have reanalyzed the data and included the stop codon region (stop codon ± 400 nt) when evaluating the distribution. Thus, we observed more than 40% m 6 A peaks located in 3'UTR and stop codon regions, and nearly 50% peaks in CDS region (Rebuttal Fig. 4, Revised Fig. 1d), which is consistent with the previous reports 8, 15 . Moreover, we had observed m 6 A enrichment around the stop codon region in scm 6 A-seq data (Original Fig. 1e, Revised Fig. 1e).
Overall, there is no significant difference in m 6 A distribution on mRNAs detected by scm 6 Aseq relative to other bulk m 6 A sequencing methods.
Rebuttal Fig. 4. The distribution of m 6 A peaks along mRNA in individual GV oocytes revealed by scm 6 A-seq. Barplot depicting the distribution of m 6 A peaks in different transcript segments. Supplementary Fig. 1e, more than 10000 m6A modified RNAs are detected in the groups of 50 and 70 GVs by bulk MeRIP-seq, while, the mounts of modified RNAs are about 5000 in the group of 25 GVs by the same way and so it is with single-cell through scm6A-seq. How to achieve the optimal balance between cell number and detection amount? Response: Thanks for the thoughtful comment. To evaluate the balance between cell number and the detection amount for m 6 A, we have merged the single GV oocytes from scm 6 A-seq into gradient cell numbers for m 6 A detection. For each gradient, we randomly selected the same number of cells for five times and merged together. The results showed that the numbers of m 6 A modified peaks (Rebuttal Fig. 5a) and RNAs (Rebuttal Fig. 5b) increase along the increasing of cell number.

As shown in
Considering the heterogeneity of single cells and the sensitivity of single cell sequencing 1 , the superposition of several cells would indeed increase the number of identified m 6 A modifications. Moreover, the sequencing depth and genome coverage might also influence the m 6 A detection in the bulk sequencing data. Therefore, we have further analyzed the ULI-MeRIP-seq data of GV oocytes with three replicates 7 to divide into gradient sequencing depth for m 6 A detection with the same method as above, and found that there was also a positive correlation between the sequencing depth and the detection amount of m 6 A (Rebuttal Fig. 5c,d).
For the results of bulk samples, it is more likely to reflect an average level of most cells, with the possibility of losing m 6 A characteristics for a small subtype of cells. Besides, the m 6 A modifications on highly expressed RNAs are more conserved and can easily be detected than those on under-expressed RNAs due to the calculation method 13 . We also obtained similar observations in this study (Original Fig 1g,h, Revised Fig. 1g,h).
Therefore, there is not a simple conclusion about the balance between the number of cells and the detection amount of m 6 A, which is related to the number of cells, genome coverage and analysis pipelines, etc. We have included this information in the discussion of the Revised manuscript (Page 13, lines 383-387). 3. Some statements lack for sufficient evidence in the article, for example, Supplementary Fig. 1h suggests that scm 6 A-seq is more sensitive than MeRIP-seq through the expression and modification of Hist4h4 gene; the integrated genomics viewer tracks of Lars2 and Mtus1 show the commonly and specifically m 6 A modified RNAs at different oocyte stages in Supplementary Fig. 4e, but it does not reflect the "far more" described in line 288. Please provide more evidence to support these descriptions. Response: Thanks for these comments. To further illustrate the sensitivity of scm 6 A-seq for m 6 A detection, we have compared the data among scm 6 A-seq and bulk MeRIP-seq 8 for GV oocytes and found that nearly 70% methylated RNAs in 2500 GVs by bulk MeRIP-seq were also identified in scm 6 A-seq (Rebuttal Fig. 6a). Moreover, we have performed the GO analysis and found that the conserved modified RNAs are annotated in cell cycle pathways, which is consistent with the previous study 8 (Rebuttal Fig. 6b). In addition, there were 50% modified RNAs only detected in single GV by scm 6 A-seq, and these RNAs are related to DNA replication and cell cycle regulation pathways, such as Hist1h4f and Mapk3 (Rebuttal Fig. 6c, Revised Supplementary Fig. 1i). To sum up, scm 6 A-seq is sensitive approach to detect both the highly conserved and low abundant m 6 A modifications which can most reflect the cellular heterogeneity (Revised Fig. 1g,h).

Rebuttal
In addition, we apologize for lacking of the statement about stage-specific modification during oocyte maturation from GV to MII. We have corrected the description in the Revised manuscript (Pages 8-9, lines 244-249), and reanalyzed and included the information about several more stagespecific modified RNAs (Rebuttal Fig. 6d-g, Revised Supplementary Fig. 4g) 4. At present, some studies have detected m 6 A modified RNAs at the developmental stages of oocyte and pre-implantation embryos, such as references 6 and 10 cited in this paper. Has the scm 6 A-seq data gained from this study been compared with these publicly available data?
Response: Thanks for this comment. As discussed in comment #2, there are much fewer reads for an individual single cell generated from scm 6 A-seq libraries than bulk samples, which results in few m 6 A modified RNAs and m 6 A peaks being detected in an individual single cell (Rebuttal Fig. 5). Thus, we have used the aggraded scm 6 A-seq data of GV oocytes to address this thoughtful consideration. Because m 6 A sequencing data in reference 6 was generated from mouse embryonic stem cells instead of oocytes 16 , we didn't include it in this comparison analysis.
To compare the scm 6 A-seq data with previous publications, we have analyzed bulk MeRIPseq data for GV oocytes with two replicates 8 , and low input ULI-MeRIP-seq data for GV oocytes with three replicates in another recent publication 7 , and compare them with the scm 6 A-seq data. We found that the modified RNAs detected by scm 6 A-seq cover nearly 70% of the modified RNAs from bulk MeRIP-seq and 55% of the modified RNAs from ULI-MeRIP-seq (Rebuttal Fig. 7a, Revised  Supplementary Fig. 1g), and moreover, the 1746 conserved methylated RNAs were enriched in the critical biological processes for GV oocytes (Rebuttal Fig. 7b, Revised Supplementary Fig.  1h). Furthermore, we found that many vital mRNAs, such as E2f1 and Bod1, are modified by m 6 A and are highly conserved among single cells (Rebuttal Fig. 7c, Revised Supplementary Fig. 1i).
Collectively, these results suggest the reliability and sensitivity of scm 6 A-seq.
Rebuttal Fig. 7 17 . This phenomenon was also illustrated in the previous publication by Kiaa1429 cKO GV oocytes and METTL3 inhibited cleavage embryos 7 . We observed that Ythdf2, which mediates degradation of m 6 A modified RNAs, expresses at a relatively higher level in the growing oocytes than the fully-grown oocytes (Revised Supplementary Fig. 3e) based on the public data 3 , which supports by the previous publication 4 about the YTHDF2 expression in the growing oocytes. Besides, we have observed a relative higher translation signal of Ythdf2 in GV than MII oocytes and zygotes (Rebuttal Fig. 8a, Revised Fig.  2e) using public data 5 . Combined with the public eCLIP-seq data 6 , we also found that the METTL3dependent m 6 A modified peaks are close to the YTHDF2 targets (Revised Fig. 2f) and more stable in Mettl3 Gdf9 cKO oocytes (Revised Fig. 2c, Revised Fig. 5d). These results suggest that YTHDF2 is the probable regulator for m 6 A mediated RNA degradation in the GV oocytes.
In another case, m 6 A modified RNAs were more stable than the unmodified RNAs during MIIto-zygote stage (Revised Fig 5b-d; Rebuttal Fig. 2; Revised Supplementary Fig. 5f). Through analyzing the published data 5 , we have interestingly observed an obvious decrease in the translation signal of Ythdf2 (Rebuttal Fig. 8a, Revised Fig. 2e), but a continued robust translation of Igf2bp2 (Rebuttal Fig. 8b, Revised Supplementary Fig. 5k) from GV to MII stage, suggesting that IGF2BP2 may play a major role in stabilizing the methylated RNAs in MII oocytes. To explain the expression change of m 6 A modified RNAs during MII to zygote transition, we have sequenced the total transcriptome of both control and Mettl3 Gdf9 cKO MII oocytes and zygote embryos. Interestingly, we found that the RNA metabolism is largely dysregulated from the Mettl3 Gdf9 cKO oocytes to zygotes (Rebuttal Fig. 8c, Revised Supplementary Fig. 5h), and a higher proportion of m 6 A modified RNAs fall in the category of hypo-stable RNAs (Rebuttal Fig. 8d, Revised  Supplementary Fig. 5i), suggesting that m 6 A contributes to the RNA instability in the Mettl3 Gdf9 cKO oocytes (Rebuttal Fig. 8e, Revised Fig. 5e). These results were consistent with the finding that m 6 A increases the stability of m 6 A modified mRNAs in normal MII oocytes during MII-tozygote transition.
Collectively, these results further illustrate that m 6 A promotes RNA degradation during oocyte growth and increases RNA stability during the MII-to-Zygote transition. More comprehensive experiments will be needed to further validate the precise role of m 6 A, such as the RNA substrate identification of reader proteins by eCLIP-seq and the metabolic analysis of m 6 A modified RNAs by RNA stability sequencing. We have included these discussions in the Revised manuscript (Page 13, lines 388-400).